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SUMMARY 


This report on the formulation of simulation equations for slung load systems is motivated by an interest 
in trajectory control for slung loads carried by two or more helicopters. A realistic, analytically simple and 
computationally efficient set of equations for use in control designs based on inversion of the nonlinear 
model is of special interest, in addition to the usual requirements for simulation and linear system analysis. 

Previous simulations and control analyses of single helicopter systems have used case-specific equa- 
tions derived from the Newton-Euler equations for rigid bodies or Langrange’s equations for general dy- 
namic systems. However, the dual lift and multilift systems are more complex, and a more efficient ap- 
proach tailored to slung load systems appeared essential to obtain tractable working equations, especially 
in the case that the suspension is assumed inelastic. 

In this report, the slung load systems are viewed as systems of rigid bodies connected by straight line 
cables or links which can be assumed elastic or inelastic. Both assumptions are of interest. The formula- 
tion for elastic suspensions has been preferred in simulations, owing to its simplicity and computational 
efficiency. The inelastic cable assumption has been preferred in trajectory control analyses, in order to 
eliminate states representing the small, higher frequency motions due to cable stretching. The conventional 
formulation for inelastic cables eliminates cable forces and cable stretching coordinates at the expense of 
increased complexity and computational requirements. These negative factors have not impeded the linear 
control analysis for single helicopter slung load systems, but are a greater obstacle for multilift system 
analysis and for nonlinear control designs. 

Equations are derived using D’Alembert’s principle and introducing generalized velocity coordinates. 
Three formulations are derived for general slung load systems. The first two generalize the previous case- 
specific results for elastic and inelastic suspensions. The third is a new formulation for inelastic suspensions 
derived from the elastic suspension equations by choosing the generalized coordinates for the elastic system 
so as to separate motion due to cable stretching from that with fixed cable lengths. The new formulation is 
significantly more efficient than the conventional one; it is computationally efficient for multibody systems 
with only a few constraints such as the slung load systems, and a compact analytical form of the nonlinear 
dynamics is obtained by appropriate selection of the generalized velocity coordinates. Further, it is readily 
integrated with the elastic suspension formulation in a single simulation, and readily applied to the complex 
dual lift and multilift systems. 

These results are applied to obtain simulation equations for the four-body dual lift system with spreader 
bar. Two other dual lift suspension arrangements of interest with three bodies are special cases and can be 
integrated with the four body system in a single simulation. The results are given as vector equations in a 
programmable form expressed entirely in terms of the natural vectors and matrices of rigid body dynamics. 
This provides a tractable set of equations for analysis and programming. 
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SYMBOLS 


A,B 

Al,L 

All, A 

Bl, ...,Bn 
Cl,..., Cm 


d 

D 

Ei(),E 2 (),E)() 

Fa 


FI, Ml 

FC1.MC1, ... 

f 

/<* 

fc 

f9 

r 

fo 

hi, hi,... 


H 

1,0 

{ia,ja,ka} 

KE 

Kl,cl, ..., Km,cm 
lol,ll, 

ml, Jl, Jn 
m 


6nx d and 6nx i r matrices of the kinematic relation v = Au + Bp. 
Column partitions of A for elastic suspension with u partitioned as 
(ul ,\) . Al, L are the Jacobians [dv/dul], [dv/d\]. 

Row partitions of A~ l in the case that u is partitioned as ( ul , \) . 

All, A are the Jacobians [dul /dv \ , [ d\ /dv] . 

Enumerated rigid bodies of a slung load system. 

Enumerated cables or links in the suspension of a slung load system. 
Number of constraints imposed by the suspension on the motion of the 
system of n rigid bodies. 

Number of degrees of freedom of the system of n rigid bodies. 

Block diagonal matrix of the system’s rigid body masses and inertias, 
diag{ml I,...,mnI,Jl,...,Jn}. 

Single axis rotational transformation about i , j , k axes, respectively, 
through angle ( ) (note 1). 

Reference frame. The subscript “a” can be any of N for Newtonian 
space; or {1 , 2 n} for rigid body axes of body Bl , ...Bn; or 
{cl , cm} for cable axes of cable Cl , ...Cm. 

Resultants of all forces and c.g. moments applied to body B\ , ..., B n . 
Resultant force and c.g. moment due to the suspension on body 
Bl,...,B n . List of resultant forces and c.g. moments on each rigid body 
(note 2). List of resultant aerodynamic and rotor forces and c.g. 
moments on each rigid body (note 2). 

List of resultant forces and c.g. moments on each rigid body due to 
the suspension (note 2). 

List of gravitational forces on the system (note 2). 

List of translational and rotational inertia reactions of the n rigid 
bodies (note 2). 

Sum of applied forces and the inertia coupling terms from /* (eq. 8). 
List of the force and c.g. moment per unit tension on each rigid body 
due to cable attachment 1 ,..M. 

List of forces and c.g. moments per unit tension due to cable 
Cl , Cm when Cl , Cm are attached to two rigid bodies. 

Also enumerated columns of the matrix H (note 2). 

Any basis of the linear vector space containing fc. 

Unit and zero matrices. Dimensions are implied by context. 
Right-handed orthonormal vectors defining a reference frame. 

Kinetic energy of a system. 

Cable spring constant and damping coefficient for cable Cl , ..., Cm. 
Unloaded and loaded cable lengths of cables Cl , ...., Cm. 

Mass and inertia matrices for the n rigid bodies of the configuration. 
Number of cables and links in the suspension. 
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M 

n 

p,-n 


Q 

r,v 


R,V 

sm 


TCI, ...TCm 
Ta,b 


(ul,X) 


X 




( )o 


O t 

()* 

eg- 

DOF 

EOM 

•,x 


Number of attachments of a cable to a rigid body of the system. 

Number of rigid bodies in the system. 

Suspension parameters which can be controlled, if any. p = (pi , ...,p*) T , 
where 7r is the number of such parameters. 

Generalized position coordinates of the system of n rigid bodies. 
q = (qi , ..., q d ) T . Also any set of position coordinates of interest. 
Configuration position and velocity, given as a list of the inertial c.g. 
position and euler attitude angles, and the inertial c.g. velocities and 
angular velocities of the n rigid bodies. The physical vectors are given 
by their coordinates in standard axes (note 2). 

Position, velocity vectors relative to inertial space. Appended numbers 
indicate specific point or line segment. 

Skew-symmetric matrix representing cross product operation with V 
for vectors referred to T a (note 1). 

Interaction force parameters for inelastic suspensions, s = (si, ..sc) T 
Cable tensions for cable Cl , ..,Cm. 

Transformation of physical vectors from frame T\, to frame F a . 

All transformations are defined from Euler angles (note 1). 

Generalized velocity coordinates of the system of rigid bodies. 
Generalized velocity coordinates u for the system with elastic 
suspension composed of 6 n — c coordinates, ul , of the system 
with inelastic suspension and c coordinates, which define 
the system motion due to suspension stretching. 

Configuration vector of inertia coupling terms in f* (note 2). 

Rigid body Euler angle triplet ($, 6, ijj) defining its attitude relative 
to inertial space, and angular velocity relative to inertial space. 

Appended numbers indicate specific body (note 1). 

Physical vector given by its coordinates in frame tF a (note 3). 

The subscipt “a” can be any of N for Newtonian space; or { 1 , 2 , . . . , n} 
for rigid body axes of bodies SI,..., Bn, or {cl , c2 , . . . , cm} for cable 
axes of cable Cl,..., Cm. 

Transpose of ( ) . 

Quantity associated with c.g. of a rigid body. 

Center of gravity. 

Degree of freedom. 

Equation of motion. 

Dot and cross product operators for physical vectors. 


Notes: 

1. Euler angles, single-axis rotation matrices, Euler angle transformations, and skew- 
symmetric matrices have the standard definitions found in aeronautical texts [1]. 
Useful formulas are listed in the appendix to this report. 
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2. The system of n rigid bodies is conveniently described by configuration vectors 
(Section 3), which are lists of vectors associated with the rigid body dynamics of the 
constituent rigid bodies. These are c.g. velocities and angular velocities, and various forces 
and moments. Standard axes are adopted for these vectors: all c.g. velocities and forces 
are given by their coordinates in inertial space, and all angular velocities and moments 

are given by their coordinates in body axes of the associated rigid body. 

3. Boldface symbols indicate physical vectors independent of any reference frame 

(R,v ,w, ....). Vectors given by their coordinates in a frame are not 
boldfaced ( R N , Vn , ....) . 



1. INTRODUCTION 


1.1 Problem Background 

Various actual and proposed slung load systems are illustrated in figure 1. Single helicopter slung 
load operations have been common since the 1950s. Such operations received extensive field use and 
development during the Vietnam war, and later research during the 1965-1975 period for the Heavy Lift 
Helicopter was focussed on the design of suspensions for stabilization of difficult loads, such as the standard 
8x8x20-ft cargo container (MILVAN). 

The use of two or more helicopters has been periodically proposed since the early success of single 
helicopter operations, using suspensions consisting of cables and spreader bars as in figure 1 ([1],[2],[3]). 
Dual lift has received limited flight test, has been used to carry payloads in a few isolated commercial 
operations, and has been advocated as an alternative to developing larger helicopters with payloads ex- 
ceeding those of current helicopters [5] or to arranging for larger capacity helicopters than those available 
in a given situation. A significant obstacle to further operational development is the complexity of system 
motion and its coordination and stabilization along any typical maneuvering flight path [6], as well as for 
precision hover [7], [8]. Progress has been hampered by lack of realistic and comprehensive equations of 
motion for use in theoretical and simulation studies of these difficulties. Tractability of the equations for 
analysis and programming, and computational efficiency, are critical factors of interest, particularly for 
control designs using recent techniques based on inversion of the simulation model [9]. While the slung 
load systems can be viewed simply as a few rigid bodies connected by cables, considerable complexity 
arises in applying the classical methods used in the previous slung load literature to the multilift systems 
when the cables are modeled as inelastic. The approach taken in this paper is to develop a systematic 
analytical formulation for general slung load systems, and to use the natural vectors and matrices of rigid 
body mechanics systematically in the applications work. The method is readily applied to the multilift 
systems and yields tractable, computationally efficient equations for the system dynamics. An alternate 
approach which avoids analysis is to consider use of one of the commercially available computer programs 
for the dynamics of general multibody systems. A third approach is to apply the previous derivation tech- 
niques to the multilift systems using symbolic digital computations to circumvent the excessive labor and 
inadequate error probabilities of extended hand analysis and programming. These alternate approaches, 
however, accommodate larger classes of dynamic systems and do not obtain the advances in tractability 
and computational efficiency achieved here by specializing to the smaller class of slung load systems. 

1.2 Equations of Motion for Slung Load Systems 

The slung load systems of figure 1 are viewed here as members of a class of systems consisting of 
rigid bodies connected by massless straight-line links which can be elastic or inelastic, and which support 
only forces along the link. These systems are characterized by the mass and inertia parameters of the rigid 
bodies and the suspension’s attachment point locations, unloaded link lengths, and link elastic parameters. 
Simulation equations are derived below for the entire class. 

The limitations of these class properties in representing the slung load systems are as follows. First, the 
bodies are assumed to be rigid; this assumption excludes from consideration elastic helicopter and rotor 
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c = # CONSTRAINTS (INELASTIC SUSPENSION) 


d = 6n-c 

= # DEGREES OF FREEDOM 

m = #OF CABLES 



SYSTEM 

n 

c 

d 

m 

1 

2 

1 

11 

1 

2 

2 

2 

10 

2 

3 

2 

3 

9 

4 

4 

2 

3 

9 

3 

5 

3 

2 

16 

2 

6 

4 

4 

20 

4 

7 

3 

2 

16 

2 

8 

5 

6 

24 

6 


Figure 1.- Examples of slung load systems. 






modes, and non-rigid loads such as half-filled tanks of liquid and long, very flexible loads. Second, cable 
mass and aerodynamics are neglected, and cable stretching is neglected in the inelastic cable model case. 
These class properties are expected to suffice for realistic representation of the low bandwidth regime of 
interest in trajectory control studies. 

Previous derivations of the equations of motion (EOMs) for single helicopter slung load analysis and 
simulation have assumed either inelastic cables ([10], [11], [12], [13]), or elastic cables ([14],[15],[16]), or 
accommodated both cable models, [17]. Previous work on the EOMs for the multilift systems has yielded 
limited results. In [18] a general formulation for these systems is given, where the bodies are assumed to 
be point masses and the cables inelastic, and [19] contains equations for the same approximation of the 
three-body dual lift system 7 of figure 1. In [7], [8] linearized equations for hover are given for the four- 
body dual lift system 6 of figure 1, assuming inelastic cables and a point mass load. Most of the previous 
work is specialized to particular suspension geometries. Exceptions are [14], [16] which account for the 
set of elastic suspensions in which all cables connect two bodies. 

The equations have been derived in the above-cited works from Lagrange’s equation for general dy- 
namic systems or from the Newton-Euler equations for rigid body motion. For elastic cables, the Newton- 
Euler equations are readily applied to all suspensions in which all cables can be considered to connect 
two bodies; this includes the multilift systems. However, in practice, cables are relatively stiff so that the 
rigid body motion with inelastic cables differs very little from that with elastic cables, and a reduced order, 
inelastic cable model is of interest in trajectory control analysis. EOMs for inelastic cables have usually 
been derived from Lagrange’s equation since the nonworking constraint forces are systematically elimi- 
nated; equations are given for a minimal set of generalized position coordinates and usually contain lengthy 
quadratic dynamic terms. In view of these factors, nearly all control analyses utilize linearized EOMs for 
inelastic cables, while simulations have been implemented for both elastic cables (e.g., [14], [15]) and in- 
elastic cables with linearized dynamics for the load-suspension motion [13]. 

The equations for inelastic cables differ considerably from those of elastic cables in form, and in an- 
alytical and computational requirements. This difference becomes acute for dual lift. The elimination of 
constraint forces in the inelastic cable derivation requires the inversion of a dx d system matrix for which 
an analytical inverse is unknown (d is the number of degrees of freedom; d = 20 for the dual lift system 6 
of fig. 1). The inversion requires an analysis of the matrix conditioning over the flight envelope of interest 
and the identification of appropriate inversion algorithms. Together with the evaluation of the large num- 
ber of quadratic dynamic terms, much more computation is required than for the elastic cable model. In 
addition, physical insight is obstructed by this matrix in the inelastic cable model, whereas the elastic ca- 
ble representation can be fully expanded in terms of the familiar vectors and matrices of three-dimensional 
rigid body mechanics. 

In the present work, simulation equations are derived by applying the Newton-Euler equations to 
each rigid body, defining configuration vectors, introducing generalized velocity coordinates, and apply- 
ing D’Alembert’s principle to the configuration of rigid bodies. This approach is similar to that taken in 
[20], [21], [22] for general multibody systems. General equations are obtained for the class of slung load 
systems, including elastic or inelastic cables, and any set of generalized coordinates. Three formulations 
are derived. The first two generalize the previous case-specific results in the slung load literature for elas- 
tic and inelastic cables. The third is a new formulation for inelastic cables obtained by appropriate choice 
of the generalized coordinates in the elastic cable formulation; that is, they are chosen to partition into d 
coordinates representing the configuration motion with inelastic cables, and c coordinates which define the 
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configuration motion due to cable stretching, where c is the number of constraints imposed by inelastic 
cables. This usually can be done by including appropriate cable velocity coordinates in the generalized 
coordinates. The result for inelastic cables is that the constraint forces are explicitly determined and ap- 
pear as an additive term in the equations of motion as in the elastic cable equations. This form requires the 
inversion of a much smaller c x c matrix (c = 4 for system 6 of fig. 1), which occurs only in the additive 
constraint force term. 

The third formulation has reduced computational penalty relative to the elastic cable equations, can be 
integrated with the elastic cable equations in one simulation, and can be expanded nearly completely in 
terms of natural vectors and matrices. It is readily applied to the multilift systems, and equations are given 
for the dual lift system 6 of figure 1. The three-body dual lift systems 5 and 7 of figure 1 are subsystems of 
the 4-body system 6, and these are readily given as specializations of the results for system 6 and integrated 
into the same simulation. The derivation and results are sufficiently brief for hand analysis and computer 
programming to be practical. 

In these results, the equations are given in terms of natural vectors and matrices with reference frames 
indicated. The generalized coordinates are comprised of natural vectors and this can be done in most ap- 
plications owing to the method of selecting these coordinates to represent an unconstrained system (elastic 
cables) with a constrained subsystem (inelastic cables). The dynamic terms can be represented entirely as 
coordinate transformations, and cross products arising from coriolis and centripetal effects and moments. 
This provides a programmable form of acceptable brevity and easy physical insight. It permits analysis and 
programming in terms of these three-dimensional vectors and matrices, and allows the devices of efficient 
coding to be applied to the vector-mechanical structure of the equations prior to treating the mass of scalar 
terms. Further expansion to scalars, if desired, is not obstructed. 

1.3 Equations of Motion for Multibody Systems 

A large body of literature on the dynamics of multibody systems has accumulated since the early 1960s 
in response to requirements in the design of machines, spacecraft, robotic arms, human motion models, 
etc., and the relevance of its formulations and results to multilift systems is of interest. The principal aim 
in this literature has been the development of general purpose computer programs to provide equations of 
motion from a minimal amount of user input data to define the multibody system. This is motivated by the 
impracticality or excessive labor of hand derivation in most working circumstances in these applications. 
Theory and analysis are given, for example, in [22], [23], [24], and surveys of computer programs which 
are equations of motion for a general system or which generate and compile symbolic case-specific code 
from user inputs are found in [25], [26]. Slung load systems differ from these applications. First, slung load 
systems with inelastic cables have only a few constraints compared to the number of degrees of freedom 
(DOFs) and can equally well be represented as unconstrained. The applications cited above are all highly 
constrained with only a few degrees of the freedom; e.g., spacecraft and robotic arms are commonly rep- 
resented as n rigid bodies with fixed orbit or base connected by n — 1 joints which permit one degree of 
freedom of relative motion, hence there are n + 2 DOFs. Consequently, a formulation in which a dx d 
matrix is formed and inverted is computationally more efficient for these applications compared to one con- 
taining a c x c matrix, while the converse is true for slung loads. Secondly, the interbody connections are 
modeled as joints in nearly all of this literature and often the generalized coordinates are predefined based 
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on the joint model. These features are advantageous for the applications cited above but their adaptation 
to cable-connected bodies is unclear. 

One code-generating program, NEWEUL, and its underlying formulation given in [20],[21],[22] avoids 
these specializations and can be applied to slung loads. The equations from this program are in the 
conventional form requiring the inverse of a d x d matrix, and can be used where this formulation is 
satisfactory. An alternate computer-based approach is to utilize symbolic computations to carry out rou- 
tine analytical steps, such as the derivatives in Lagrange’s equation and routine code optimization. The 
general possibilities of applying MACSYMA [27] for this are discussed in [28] and the approach is used in 
[8] to derive linearized hover equations for the dual lift system 6 of figure 1 from Lagrange’s equation. In 
general, these computer-generated codes are fully expanded to scalar components and the resulting pages of 
expressions are unsuited to hand analysis or physical insight in larger problems, but substantial advantages 
are obtained where this is not a consideration. 

For slung loads, the present paper provides a method which makes hand derivation, analysis, and pro- 
gramming feasible for the dual lift and multilift systems, and a new form of the equations which improves 
computational efficiency of the inelastic suspension case over previous forms. Results are given for the 
dual lift systems. A further technical paper is in progress which expands the present material to include 
EOMs for additional slung load systems in figure 1, EOMs for degenerate body (point mass, rigid rod) 
approximations, and linearized EOMs both for general slung load systems at general reference flight con- 
ditions and for the dual lift system with spreader bar. 

2. SYSTEM DESCRIPTION 

The systems of interest consist of one or more helicopters which support a load (more than one load in 
some instances) by means of a suspension. The load to be supported by the helicopters due to gravity, load 
acceleration, and aerodynamics is dominated by gravity for typical slung loads and nominal trajectories. 
The suspension can include spreader bars of non-negligible weight (e.g., the multilift systems) and other- 
wise consists of straight line links made up of cables, usually of nylon webbing, along with hooks, rings, 
bars, isolator springs, or other hardware [29]. Suspension designs with controllable geometry obtained by 
active cable winching and attachment point movement have been proposed for load attitude stabilization 
in [30], [31], and are included in the present formulation. 

These systems contain n rigid bodies, B 1 , .. .Bn (helicopters, load, spreader bars). The remainder of 
such a system (termed the suspension hereafter) consists of m straight-line links which support only force 
in the direction of the link (only tension in the case of cables), and has negligible mass and aerodynamic 
force compared to those of interest in trajectory control. Every link is connected at both ends to a rigid 
body or to another link. The links can be modeled as inelastic, in which case c<m holonomic (position) 
constraints are imposed on the motion of the rigid bodies, and the system has d = 6 n — c degrees of 
freedom. Values of c, d are listed in figure 1 for the systems shown there. For example, in system 3 of 
figure 1 the load and 4 inelastic cables form a single rigid body with a point held fixed to the helicopter, 
so that c = 3 . In classical mechanics, the motion of such a system is usually represented by d EOMs for d 
independent accelerations with the nonworking constraint forces systematically eliminated. Additionally, 
equations for a minimal (d) or nonminimal set of motion variables can be given with the constraint forces 
present, but expressed in terms of the applied forces and states. Alternatively, the links can be modeled as 
elastic due to cable or isolator spring stretching, in which case there are 6 n DOFs and the corresponding 
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6 n EOMs require calculation of the forces applied by the suspension to the bodies from the suspension 
stretching and cable elasticity parameters. 

An examination of cable and suspension elasticity and its effect on rigid body motion is found in [32]. 
Cable stretching under tension is usually modeled as that of an undamped spring with damping supplied 
by the aerodynamic resistence of the attached bodies. Cables tend to be stiff, but the suspension design 
must avoid an upper bound given by resonance with the helicopter rotor frequency (around 4-5 hz), where 
a divergent pilot-induced oscillation vertical bounce mode has been observed near hover. The net result 
is natural frequencies of practical suspension links around 2.0 - 2.5 hz. This frequency is sufficiently high 
to be disjoint from the frequency range of interest in trajectory control (around 0.25 hz). The correspond- 
ing mode is one of rapid and significant cable tension variations, but with small stretching excursions so 
that the rigid body coordinates are almost unaffected [16]. Both elastic and inelastic suspension models 
are of interest in trajectory control. For practical suspensions, simulations can employ the simpler, more 
computationally efficient nonlinear equations of the elastic model. If suspension stiffness were signifi- 
cantly greater, then difficulties of numerical stability and ill-conditioning would arise in real-time digital 
simulation of the higher frequency, lower amplitude cable stretching motion. In control designs, practical 
suspensions are approximated as inelastic to eliminate feedback of states with negligible influence on the 
rigid body motion. If suspension stiffness were significantly lower, then the lower frequency, higher am- 
plitude motions due to cable stretching would be of interest in trajectory control. The present report studies 
simulation equations for both elastic and inelastic suspensions. 

The simulation of cable collapse is an application detail outside the present scope. If a cable collapses, 
the resulting system is still a member of the class and can be simulated. However, practical suspensions 
are designed and operated such that cable collapse does not occur except during large unstable excursions 
from the nominal configuration, or a redundant cable can collapse under ordinary motions but leaves the 
constraints imposed by the suspension invariant. 

3. EQUATIONS OF MOTION OF GENERAL SLUNG 
LOAD SYSTEMS 

3.1 Configuration Vectors 

Physical vectors are referred to inertial or body axes reference frames in the following. Translational 
motion and forces are given in inertial coordinates, and rotational motion and moments in body axes. 
The reference frame is indicated by a subscript, which is N for inertial space and i E { 1 , 2 n} for 
body axes of the ith body. Body axes components of translational velocity and motion variables relative 
to a reference body are commonly used in slung load simulations, and are readily introduced later when 
selecting generalized coordinates for an application. 

It is convenient to use configuration vectors which define the motion and forces of the » rigid bodies, 
whose masses, inertias, c.g. translational motion, Euler attitude angle triplets, and angular velocities rela- 
tive to inertial space are denoted by (ml , J1 ,RI*, VI *,ck1,u> 1 (ran, Jn, Rn\ Vn*,atn, wn). The 
configuration vectors of position, r, and velocity, v, are defined as lists of the rigid body c.g. positions and 
Euler attitude angles, and the rigid body c.g. translational and angular velocities: 
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/m* N \ 


(VI *N\ 

Rtf n 
al 

v = 

Vn' N 

wli 

\ an / 


\ UTln ) 


( 1 ) 


Let / be a list of the resultant forces and moments applied to each rigid body and let fg, fa, fc be 
corresponding lists of the gravitational forces, the aerodynamic and rotor forces, and the constraint forces 
due to the cables, respectively: 


f = fg + fa + fc 


/ml g N \ 


FA1 n \ 


FC 1 * \ 

mng N 

f a - 

FAn N 

fc = 

FCn N 

0 

MAh 

MCh 

0 ) 


\ MArin ) 


\ MCrin / 


where FAi , MAi are the sum of the aerodynamic and rotor forces applied to the ith body and the sum of 
their moments about its c.g.; and FCi , MCi are similar force and moment sums due to all cables acting on 
Bi. The applied forces, f, depend on r, v, and the helicopter controls; any dependence of the aerodynamics 
on the acceleration, v, is assumed negligible. 

Last, let f* be a list of the inertia reactions of the n bodies: 


/ ml VIn \ 


/ 


0 


\ 


mn Vun 
J 1 wli 


0 

S(uh)Jl uh 


= -D v - X 


(3) 


\ Jn UUn ) 


\ Siwrin) Jn / 


f* is abbreviated to the form f* - —D v — X where D is block-diagonal with masses and inertia matrices 
along the diagonal, v is the configuration acceleration, and X contains coriolis terms due to the use of body 
axes components of rotational motion. 


3.2 Kinematics 

The systems of interest consist of n rigid bodies connected by m cables which impose c constraints on 
the motion of the rigid bodies; they have d = 6 n — c DOFs, where c = 0 for elastic cables. These are 
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holonomic systems; that is, the constraints imposed by an inelastic suspension can be given as functions of 
position only. These constraints are usually time-invariant, but in the special case of active cable winching 
or attachment point movement these have explicit dependence on time. To accommodate this case, the 
controllable geometric parameters of such a suspension, p = ( pi , . . . , p „) T , are included in the development 
below. 

For holonomic systems with d DOFs there exist d generalized position coordinates, 
q = ( Q\ i •• •> Qd) T , which suffice to locate all points in the system and also the configuration position: 

r = r(q,p) (4) 

and d generalized velocity coordinates which suffice to define all inertial velocities of the system: 

u = U(q,p)q (5) 

The configuration velocity is related to u by a linear expression of the form: 

v = A(q,p)u+ B(q,p)p (6) 

Here, U is a nonsingular dxd matrix; it can be unity, but velocity coordinates different than q are commonly 
useful in applications. Note that v is asserted in equation (6) to be linear in u,p. This follows from the 
usual linear relationship v(f) from rigid body kinematics (appendix) 

V i* N — Ri* N 1 . , 

• Txss \ ■ - 

uiy - W{ai) at J 

and equations (4) and (5). A is a 6 nx d matrix. For inelastic suspensions, A expresses the constraints 
by confining that part of the configuration velocity due to u to the d-dimensional linear vector subspace 
defined by the columns of A. If p = 0 then this subspace is tangent to the configuration trajectory, r(t) . 

3.3 Equations of Motion from D’Alembert’s Principle 

D’Alembert’s principle applied to a rigid body, Bi, of a holonomic system yields the following equations 
for the inertial motion (e.g., see [33]): 

( FiN — mi Vi*N) T dRi*N = 0 

[Mi, — Ji tin',- — S(uii) Ji = 0 

where Fi , Mi are the total applied forces and moments about the c.g. of body Bi, and dRi *, d ai are the 
c.g. position and rigid body attitude changes due to any perturbation of the generalized coordinates, dq, or 
equivalently, due to any small motion, udt : 

dRi* N = [^~~\udt 
au 

dwu 

daii = [^]udt 
du 
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The velocity Jacobians [dVi* N /du] and [duii/du] are the appropriate rows of A in equation (6) and the 
result is valid for arbitrary u. The same result for all n rigid bodies of a slung load system can be combined 
to get the following d equations of motion: 


A T (f+n= 0 ( 7 ) 

This result can be obtained almost as easily from Kane’s equations [33], since A T f and A T f* can be 
recognized as lists of the generalized active and inertia forces for each DOF. Lagrange’s equations have 
been used in much of the previous literature on slung load systems with inelastic suspensions, and can 
be applied to the general slung load system where kinetic energy is the sum of the kinetic energies of the 
constituent bodies: 

KE = \v t D v = i(<f [U t A t DAU] q + 2q T [U T A T DB]p + p T B T DBp ) 

For holonomic systems, A T f and A T f* are equal to lists of the generalized applied and conservative 
forces, and the derivatives of the kinetic energy [33], respectively, in Lagrange’s equations. D’Alembert’s 
principle is used here for slung loads in order to maintain contact with the underlying rigid body dynamics. 
One consequence is the compact, well-structured formulation of the second-order velocity dynamics in 
terms of the rigid body velocities, v, seen in X (eq. 3) and in the remaining such terms in the dual lift 
example given later. Lagrange’s equations lose contact with v in KE and generate the second-order velocity 
dynamics in terms of the generalized velocity coordinates, q, from derivatives of KE. The result is a large 
number of scalar terms with no apparent structure, and the number of such terms grows rapidly with the 
number of rigid bodies. 

To obtain the simulation equations, differentiate equation (6) with respect to time, introduce this in f* 
(eq. 3), and solve equation (7) for u: 

ii = [A T DA]~ 1 A T [fo+ fc] (8) 

where 

fo = fg + fa - DAu - X- D(Bp+Bp ) 

The configuration vector fo denotes the combined external forces, second-order velocity effects due to the 
choice of coordinates u, i>, and the inertia reaction of the configuration to p( t ) . 

If the cables are elastic, then A is a nonsingular 6 n x 6 n matrix and: 

u = A~ l D~ l [fo + fc] (9) 

In the case that p = 0 and we choose u = v, then A - I and the result is identical to the Newton-Euler 
equations applied to each body: 

i> = D- l [fg + fo-X + fc] 

This result can be given by inspection of the system and has been the starting point in most work which as- 
sumes elastic cables. Equation (9) generalizes this to allow any choice of generalized velocity coordinates; 
for example the use of cable velocity coordinates in u provides a convenient and well-conditioned calcu- 
lation of cable lengths and directions for use in the calculation of the cable forces, f c . It is remarked here 
that the matrix A~ x in equation (9) represents the kinematic relation u(v) and can be given analytically 
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from the kinematics as readily as the matrix A representing the reverse relation v( u) . It is unnecessary to 
perform matrix inversion to obtain A ~ x . 

The constraint force, /c, in equation (9) can be given as a sum of forces and moments applied by the 
suspension at each attachment point: 

fc^hjTCj 

;=i 

where j enumerates every attachment of a cable to a rigid body, hj is a configuration vector defined in the 
next section, and TCj is the cable tension which is given by the spring model of the cable as 

TCj = moi{0 , Kj(tj - Zoj) + cj ioj) j - 1,2 M 

where {ioj, Kj , cj } are the unloaded cable length, and cable spring and damping constants. Cable damp- 
ing, cj j- 0 , is introduced in [16] but otherwise this has been neglected in simulations with elastic cable 
models. The spring constants if, care related to the natural mode parameters by w* = K/m,2£uj n = c/m, 
where m is the load mass supported by the cable. In simulations with elastic cables, initial values of cable 
tension and stretch can be calculated in cases with nonredundant suspensions from the solution for the 
constraint forces, / c , for inelastic cables given below. Alternatively, the configuration can be allowed to 
settle from an approximate initial arrangement, possibly with the aid of cable damping as a settling device. 

3.4 Equations of Motion for Inelastic Suspension: Method 1 

If the cables are inelastic then the constraint forces, fc, drop out of equation (8). This result is shown 
as follows. First, enumerate the cable attachments 1,2,. ..,M at attachment points RI , R2 , ...,RM on their 
corresponding bodies Bi( 1) , Bi( 2) , . ..,Bi( M). One or more cables are attached at an attachment point 
and every such attachment of a cable to a body is numbered. The constraint force on the configuration is 
then (see fig. 2): 

/c = X>;rC; (10) 

;=i 

where TCj is the cable tension. The j th cable attachment at point Rj on body i(j) applies a force and 
moment to Bi(j) given by: 

FCij = kcj TCj 
MCij = (Ri’j x kcj) TCj 

where kcj is the cable direction outward from the body and Ri *j is the attachment point to c.g. moment 
arm (Rj - Ri*). Thus, hj in equation (10) is a configuration vector whose nonzero elements are kcj and 
( Ri *j x kcj ) , corresponding to the constraint force and moment on Bi(j) due to the j tk cable attachment, 
per unit tension. 

Second, enumerate the cables and links which constitute the suspension, Cl, C2,...,Cm. Each end of 
a cable or link is attached to either a rigid body or to another link. In the special case that all cables are 
attached at both ends to rigid bodies then all terms in equation (10) can be combined in pairs, with each 
pair corresponding to a single cable: 


i 
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Figure 2.- Suspension forces. 


/c = £ Hj TCj (11) 

;=i 

where the nonzero elements of Hj are kcj , -kcj , ( Rj — Ri *) x kcj , — ( R j ' — Ri '*) x kcj corre- 
sponding to the forces and moments on bodies Bi, Bi' to which the cable C j is attached at points Rj , Rj ' 
(see fig. 2). Here, kcj is the direction of the line segment ( Rj ' — Rj ) . From the definitions of Hj , v it fol- 
lows that v T Hj combines the rigid body c.g. and angular velocities to give the difference in cable endpoint 
velocities along the cable direction: 


v T Hj = Vi* • kcj — Vi * • kcj + wi • (Ri*j x kcj) — wi' • (Ri *j' x kcj) = ( Vj — Vj') • kcj = —tj 


that is, it is the cable stretch rate. Recalling equation (6), then u T A T Hj is the cable stretch rate due to u. 
Since all cable lengths are invariant for arbitrary u in an inelastic system, then: 


A T Hj = 0 ; = 1,2,. ...m 


and 

A T fc = 0 

Equivalently, the constraint forces do no virtual work, f(F8r = 0 . All systems in figure 1 can be arranged 
such that every cable connects two rigid bodies. 

For more general suspensions, fc is given by equation (10) and v T fc is a linear combination of 
the attachment point velocities. To this combination, first add and subtract the velocities of the cable 
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interconnection points in the cable directions, and then apply the force balance condition to the linear 
combinations of cable forces at these interconnections. The result is: 

Mm to 

t> T /c = £Vj • kcj TCj = £( Vj - Vj') • kcj TCj = -££; TCj (12) 

;=1 ;=1 J =\ 

where the second and third sums are taken over all cables and {lj } are cable lengths. Consequently, if the 
suspension is inelastic then all cable lengths are invariant for arbitrary u and: 


A T fc = 0 


(13) 


Equation (8) for inelastic suspensions is now: 

u = [A T DA]~ l A T fo (14) 

Here, [A T DA] is a positive definite d x d matrix, D contains the system mass-inertia parameters, A 
contains the system geometry and constraints, and fo (defined by eq. 8) contains the term DA u + X, 
which is second-order in velocity coordinates, and the inertia reactions to p. It is unnecessary to calculate 
cable tensions in this result. 

The principal difficulty in applying equation (14) is the required inversion of [A t Dj 4]. A general 
analytical inverse is unknown. For slung loads, d is relatively large (near 6 n) and this inversion is com- 
putationally costly compared to the elastic cable case (eq. 9), which requires no inversion. Some analysis 
of its numerical conditioning over the flight envelope is needed before inversion can be confidently im- 
plemented. Additionally, the matrix inversion obstructs further expansion and rationalization in terms of 
three-dimensional dynamics while equation (9) can be fully expanded in applications in terms of natural 
vectors and matrices to obtain the sensitivities of the elements of u to the forces and moments. A Gauss- 
Jordan elimination procedure to reduce the matrix to a block diagonal form was considered in initial efforts 
to reduce the inversion to the inverses of smaller blocks given in terms of three-dimensional vectors. The 
procedure was applied in each case after defining the coordinates, u, and identifying some subdivision 
of u into natural vectors or groups of coordinates with which to associate the blocks of the diagonalized 
form. For typical choices of u, these blocks usually represented the rigid body inertia matrices modified 
by the constraint moments. Analytical results were obtained for systems 1, 2, 3, 5 of figure 1, but the labor 
required increased rapidly with the number of DOFs and the procedure appeared impractical for the dual 
lift system 6 by hand. No generalization of the “effective inertia” representation of the blocks or of the 
simplification required at each step of the reduction was found, so the procedure was ad hoc and specific 
to each system and choice of coordinates, and could not be given effectively as a MACS YMA algorithm 
for analytical reduction. 

A second form of the simulation equations requiring the inverse of a much smaller cx c matrix, was 
obtained after restricting the generalized coordinates appropriately for slung load systems and is described 
next. The result can be applied with much less labor and likelihood of error than the procedure considered 
above. 

Equation 14 or analogous forms are implemented in many of the general multibody programs where 
d is relatively small in typical applications. Numerical inversion of positive definite matrices can be car- 
ried out by various algorithms based on the Cholesky (square root) decomposition ([34]) previously ap- 
plied in the field of estimation. These are computationally efficient and resist numerical instabilities. The 
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conditioning of the coefficient matrix [ A T DA ] depends on A in equation (14) or, equivalently, on the 
choice of coordinates, u. In the multibody programs these coordinates are often preselected based on the 
joint model of interbody connections, and these appropriate coordinates tend to result in well-conditioned 
coefficient matrices. In the present work, the interbody connections are suspensions comprised of cables; 
it is left to the applications phase to determine in each case what constraints are imposed and what choice 
of coordinates best represents the constrained system motion. 

3.5 Equations of Motion for Inelastic Suspension: Method 2 

The approach in this section is to select the generalized velocity coordinates for the elastic system in 
equation 9 so as to separate motion due to cable stretching from the remaining motion with invariant cable 
lengths, and then study the results when the stretching motion is nulled. A solution for the constraint force 
and new simulation equations for the inelastic system are obtained. 

First, let the generalized velocity coordinates of the elastic system, u, be comprised of 6n — c coor- 
dinates, ul , of the system with invariant cable lengths and c length rates, A, which suffice to define the 
motion due to cable stretching. In general, if the c independent position constraints imposed by an inelastic 
suspension are given by (Al(r) = 0, Ac(r) =0} then A can be taken as ( A 1 , Ac) r . For slung load 
systems, A can usually be taken as cable lengths and the complete set of coordinates (ul , A) taken as the 
c.g. velocity of a reference body, the angular velocities of all rigid bodies, and the cable angular velocities 
and stretching rates. 

Next, substitute the partitioned u in equations (6) and (9): 


v = Au + £p=Al ul + L A + £ p 

(15) 


(16) 


where fo is defined in equation (8); A1 , L are the 6 n— c and c columns of A which, respectively, define the 
contributions of ul , A to v; and AI1 T , A r are the 6 n- c and c rows of A~ l which respectively define the 
relations ul( v) and A( v) . As noted earlier, A -1 can be obtained without matrix inversion since it defines 
the relation u( v) for the elastic system and can be given from the kinematics as readily as A, which defines 
the relation v(u) . 

Equation (16) gives the simulation equation for systems with elastic suspensions in terms of the coordi- 
nates ( u 1 , A) , where u 1 leaves the cable lengths invariant. The influence of cable stretching motion on u 1 
can be viewed by entering the derivative of the partitioned generalized velocity coordinates equation (15) 
in /* and rederiving equation (8). The first 6n — c equations can be arranged as 

ul = [Al T D Al]" 1 A\ t [fo - D(L\ + LA)] 

where fo is as defined with equation (8) except that Al ul replaces A u. As in equation (14), fc drops 
out ( A 1 T fc = 0 in view of equation 1 3) and the result differs from the inelastic system equations (eq. 14) 
only in the presence of the configuration acceleration due to elastic stretching. The effect of the stretching 
motion on ul depends on the cable spring constants: for a fixed disturbance, the extremes of A decline 
with increasing spring stiffness, while the extremes of the term in A remain fixed in magnitude, but A(t) 
increases in frequency [16]. 
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For inelastic suspensions, X = 0 and equation (16) gives c equations from which to calculate the 
constraint forces: 

0 = A T D~ l (fo+fc) (17) 

The constraint force can be given in terms of c independent parameters; that is, it can be given in the form: 


(«) 


fc=[Hl,...Hc] 


= H a 


W 


where s is arbitrary, {#,} are configuration vectors, and rank [ H] = c. From equation (13), A\ T H = 0 , 
and from the construction of A ~ l , A\ T A = 0 . Therefore, the columns of H and A are both bases of the 
same linear vector space, and A can be used to define the constraint force: 

fc = A s (18) 


where A T is the Jacobian, [dX/dv], and s has units of force. Introduce equation (18) in equation (17) and 
solve for a: 

a = -[A T D- l A]~ l A T D- { fo (19) 

Further, A can be replaced in equations (18, 19) by any other convenient basis, H , of the space containing 
fc. For example, in the special case that every cable connects two bodies and c = m, then s can be taken 
as the c cable tensions with the basis vectors, HI , He, as defined by equation (11) above. 

The cable tensions are related to s, and can be uniquely determined from s if the suspension is not 
redundant. However, the constraint force applied to the configuration, fc, can always be calculated. A 
suspension can be separated into disjoint sets of interconnected links. Each such set imposes one holonomic 
constraint; if the number of sets is c then all cable tensions can be uniquely determined from s, but if it 
exceeds c then the constraints are redundantly imposed. In the special case that all cables connect two rigid 
bodies, then each cable is a disjoint set and all cable tensions can be found only if c = m. 

Finally, the simulation equations are: 

u = A-'D~'(fo+ A a) (20) 


or, for the inelastic DOFs: 


til = An T D-'[I-A[A T D-'A]-'A T D-']fo (21) 

This result has several advantages over equation 14 for the slung load systems. First, the leading coefficient 
matrix, AI 1 , is known analytically. Second, equations (19), and (20) require the inversion of a c x c matrix, 
[A T D~ l A], where c = 1 for system 1 and c = 4 for the dual lift system 6, compared with the required 
inversion of an 11 x 11 and a 20 x 20 matrix, respectively, for these two systems when equation (14) is 
used. In many cases, the cable tensions can also be generated by equation (19). 

Computational efficiency in calculating the dynamics can be compared by counting the number of 
multiplies required to generate ti, given A, A~ x , D, D ~ x , /o . These are shown in table 1 for systems 1, and 
6 of figure 1. 
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Table 1. Number of multiplies to calculate ti 



System 1 

System 6 

General Formula 

Eq.(9) 

Eq. (14) 

Eqs. (10-20) 

60 

616 + g( 1 1 2 ) 

119 

264 

3430 + g(20 2 ) 
868 + g(4 2 ) 

6n(6n - i +2) 

d((6n - i)^|^- + 12n -t-d) + g(d 2 ) 
d(6n -i) + 3n(c 2 + 9c + r) + c 2 + g(c 2 ) 


Here, g() refers to the number of multiplies and divides required for the Cholesky inversion, which 
increases with the square of the matrix size ([34]). The general formulas (derivations omitted) include 
savings due to generic zeros, ones, and matrix symmetries. The number of coordinates u, v which is 
identical is represented by t in these formulas, and accounts for generic zeros and ones. In most cases, 
u can be selected such that * = 3n + 3. Additions are omitted from the operations count for brevity 
but this omission does not affect the conclusions. Equation (14) is representative of previous single-case 
formulations of the slung load dynamics with inelastic cables, as well as that in most of the multibody 
literature (e.g., [22]). As seen, equations (19) and (20) provide a significant reduction in the computational 
requirements to represent these dynamics, and the same conclusion applies to all systems of figure 1. 

3.6 Procedure for Generating Simulation Equations 

A procedure for applying and implementing the previous results is outlined in figure 3. Both equa- 
tion (16) for elastic suspensions and equations (19) and (20) for inelastic suspensions are included. 

First, (figure 3a) analysis is required to determine c, define appropriate generalized velocity coordinates 
(ul , \) , and obtain expressions for A, A~ x , A from the kinematics, and also expressions for B, B in the 
case that p ^ 0. Secondly, the simulation equations (16) or (19, 20) can be implemented as shown in 
figures 3b and 3c. For elastic cables, fc is obtained from {hj, TCj} which are defined with equations (10) 
and (11) above, and for inelastic cables, fc is obtained from equations (1 8) and (19). Subsequently, the rigid 
body kinematics v, r are obtained in the usual way or, more generally, this can be expanded to any set of 
velocity and position coordinates of interest, q, q after defining q(u,p ) from the kinematics and geometry. 
For inelastic cables, the coordinates \ and their equations can be eliminated, since \ is theoretically zero. 
Alternately, all 6 nequations can be retained to permit simulation of both cable models; for inelastic cables, 
the computed X indicates numerical accuracy. 
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c 


number of holonomic constraints on configuration for inelastic suspension. 


ul 6n-c generalized velocity coordinates for system with inelastic suspension. 

X c generalized velocity coordinates defining suspension stretching motion. 

u 6n generalized velocity coordinates for system with elastic suspension 

u=(ul r ,X r f 

A(r,p) the Jacobian [dv/du] obtained from the kinematics. A=[A1IL] 

A -1 (r,p) the Jacobian [du/dv] obtained from the kinematics. A -r =[AIllA ] 

B(r,p) the Jacobian [dv/dp] (Required only if p y 0). 

(a) Quantities required a priori. 

fa = fa(r,v,6) (6 = helicopter controls) 

fo =fg + fa-DAu-X - D(B p + Bp) 

f c hj max{0 , Kj {lj - ioj)} elastic cables 

— A[A r Z) -1 A] _1 A r Z> -1 /o inelastic cables 

sf = D" 1 (fo + fc) 

u =A _1 sf 

u = / udt 

v =Au+Bp 

f = r(v) or q- q(n,p) 

r =/r dtorq=/g dt 

(b) Simulation equations. 

Figure 3. Procedure for simulation of slung load systems 
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Figure 3(c).- Simulation flow diagram. 

4. DUAL LIFT SIMULATION EQUATIONS 

Simulation equations for the dual lift system 6 of figure 1 are derived next. Equations for systems 5 and 
7 of figure 1 are are inelastic and are part of the load rigid body, B4, and that the system has no controllable 
geometry (p = 0) . The derivation uses formulas from the appendix without further citation. 

The first task is to define generalized velocity coordinates for the system with elastic suspension which 
separate system motion with invariant cable lengths (20 coordinates) from the motion due to cable stretch- 
ing (4 coordinates). To do this, observe that the inertial c.g. position and velocity of every body in the 
system can be given from the c.g. position and velocity of a reference body (taken as Bl) plus any se- 
quence of straight line segments in the system which connect the two c.g.’s and the velocities of these 
segments. For example: 

R4* = Rl* + Rl*l + R13 — R3*3 + R3*5 — R4*5 \ (m 

V4* = VI* + wl x Rl*l + V13 — w3 x R3*3 + V3*5 - w4 x R4*5 J { } 

The segments are c.g. attachment point line segments ( Rl *1 , R3 *3 , R4 * 5 ) or cables ( R13 ) or defined 
from cables ( R3 *5 ) . The velocities of the fixed length segments ( Rl *1 , R3 *3 , R4 *5 ) are given from 
their angular velocities as seen in equation 22. Thus, it suffices to define generalized velocity coordinates 
which separate the cable velocities V13 , V24 , V3 *5 into motion due to cable stretching and rotational 
motion with invariant cable lengths. 
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Bl. HELICOPTER 1 
B2. HELICOPTER 2 
B3. SPREADER BAR 
B4. LOAD 

Cl. TETHER CABLE 1 
C2. TETHER CABLE 2 
C3. BRIDLE CABLE 3 
C4. BRIDLE CABLE 4 


ml, Jl, Rl*lj 
m2, J2, R2*2 2 
m3, J3, R3*3 3 , R3*4 3 
m4, J4, R4*5 4 

tl 

12 

13 

14 = 13 


Figure 4.- Dual lift with spreader bar. 
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The tether cable velocities are treated as given in the equation table, table 2a. A cable axis, T c , is 
defined for each cable with its A: -axis along the cable, and oriented inertially by pitch and roll angles from 
which the transformation to cable axes, T C)N , is given. These angles are small when the cables are near the 
vertical. Expressions for the cable-axes components of the cable velocities VT3 c i, V24 c2 are given from 
derivatives of the appropriate position coordinates as shown, and are seen to be appropriate generalized 
velocity coordinates which separate rotational and stretching motion. Linear velocity components are used 
in preference to cable angle rates for simplicity in the relation u( u) given later. 

The velocity contributions of the bridle cables, C 3 , C4 are formulated in table 2b. These cables form a 
triangle with the spreader bar longitudinal axis, i3 . Axes, are attached to the plane of the triangle. The 
triangle plane can have any roll angle, <f> t , about i3, and R3*5 has components (£3x,£3z) in this plane. 
To isolate the velocity coordinates ( <£ t , 13 x, 13 z) in an expression for V3 *5 define a reference frame, Fy, 
from the spreader bar Euler angles ( # 3 , ^ 3 ) and then evaluate the inertial velocity of R3 *5 as the sum of its 
velocity relative to Fy, denoted as VT , and the effect of the inertial rotation of Fy, denoted as A VT . The 
three components of VT t are the desired generalized velocity coordinates. A VT depends on the angular 
velocity of Fy, which can be given from the spreader bar angular velocity as shown in table 2b. 

Finally, the generalized coordinates, u, are assembled and related to the configuration velocity, v, by the 
matrix, A, given in table 2c. For simplicity of A, u contains as many of the rigid body velocities as possible, 
in addition to the required cable velocity coordinates. Here, u includes the c.g. velocity of a reference 
body, V 1 *n , and the angular velocities of all bodies. This can be done for most systems of figure 1 and 
the corresponding rows of A are from the unit matrix. The remaining rows generate V2*n,V3*n,V4*n 
from kinematic expressions analogous to equation ( 22 ) and contain only coordinate transformations and 
skew-symmetric matrices representing cross products with the c.g.-attachment point line segments in the 
coriolis velocity terms. 

The rows of A~ x which are not from the unit matrix generate the cable velocities, V13 c i , V24^,VTt, 
as differences of the attachment point velocities obtained from the rigid body velocities of the two bodies 
connected by the cable, and contain the same matrix objects found in the submatrices of A. 

The configuration forces are assembled next. First, the applied forces and inertia coupling terms, fo, are 
assembled in table 2d as the sum of applied forces due to weight (ml g , . . . , m4 g ) , the forces and moments 
due to the aerodynamics and helicopter rotors ( FA1 , . . . , MA4 ) , and the inertia coupling term, X + DA u. 
The submatrices of A are defined in table 2c and their time derivatives are obtained principally from the 
derivatives of various coordinate transformations. A general formula for transformation rates is given in 
table 2d along with expressions for the cable axes angular velocities and the results for the terms from A u. 
These terms are seen to be principally coriolis and centrifugal accelerations and are compactlyexpressed 
in terms of the rigid body velocity coordinates. The time derivative of BT is obtained routinely, but it is 
omited from the equation summary for brevity. The vector elements of fo are denoted F01 , ..., M04 for 
brevity in the later equations. 
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Table 2. Simulation equations for a dual lift system with spreader bar 


(a) Velocity coordinates for tether cables ((71,(72). 

T c = {ic,jc, kc} c 6 {cl, c2} 

Tc,n = E 1 (<f> c ) E 2 (9 c ) 

kc = costfic sin0 c iN — ain4> c jN + coa<j> c cos0 e kN 
RV$n — £\ kclff 

R24m = {-2 kc2pi 

V13 d = T c i t N R13 n = (£i 0d cos<f>d, —£\ 4>c\ , £\) T 
V24 c2 = T C 2,N R24n = (£ 2 #c2 COS(f> c2) —l 2 rf>c2 > h) T 

(b) Velocity coordinates for bridle cables (C3,C4). 

Tt - {i3,jt,kt} 

i3 =spreader bar longitudinal axis 
kt = perpendicular to i3 in plane of (73, (74 
R3*b t = (f3*,0,f3z) T 
Tv ={i3,j3',k3'} 

Tv,n = E 2 (0 3 )E 3 (^ 3 ) 

Tt, v = Ei{tf> t ) 

V3*5 1 = T tlN R3*5 N = VT t + A VT t 
jR3*5 3 i = r t)3 . i23*5 t 

VT t = T t ,v R3*bv = {£3x, -£3z j> u h z) T 

AVT = — R3*5 x w3' = -R3*5 x |j3' 0 3 + kN ^ 3 ] = -R3*5 x |j3' j3' + kN k3'/cos0 3 ] • w3 

( 0 £3z cosA.<f> £3z sinA<f> \ 

0 — £3x sinA4> + £ 3z sin<f> 3 tan0 3 £3x cosA<f> + £3 z cos<f> 3 tan0 3 I o*3 3 
0 — £3x cosA(f> —£3x sinA<j> J 

where A<j> = <j) t — <f>z 
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Table 2. (continued) 

(c) Configuration kinematics. 




An = -T Ntl S(Rl'U) 

An = T n .iS(R2"2 2 ) 

4 37 = Tn i3 S(R3*3z) 

A' 27 = T n< 3 S(i23*4 3 ) 

■4.27 = -4-37 — 4.27 = —Tn,3 S(R34 3 ) 

4 47 = 4 37 + Tx <t BT 

4 4 8 = T;v, 4 5'(i24*5 4 ) 


#25 = -T c \,nA-25 
#36 = T c 2,nA26 
B 2 7 — — 2c1,n4 3 7 
-®37 = ~T C 2,N A 27 
B\ 8 = — T t) AT 4 4 8 
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Table 2. (continued) 

(d) Applied forces and inertia couplings. 


f° = fg + fa — X — D Au = 


(F 01 N \ 
F 02// 

fqz n 

F 04^ 
M Oli 
M 02 2 
M03 3 
\M04 4 / 


/ 


W] 5AT + .FAljy 

9N + FA2n - ra 2 [Tjv,ci V13 c i - Tn,c 2 1^2402 + A 2 s o;li 

+-4-26 0>2 2 + A2 7 0>3 3 ] 

m 3 sw + FA3n — m 3 [f n,ci V13 c i — A25 o,l 3 + A 37 o;3 3 ] 

m 4 9N + FA4n — ui4[2V,ci V13 c i — VTt + A 2 s o;li + A 47 o;3 3 + A 4 g o>4 4 ] 

A/ Ali — 5(wli) J1 a>li 

Af^4.22 — *S(u?22) «/2 u>22 


\ 


I MA3 3 -5(o/3 3 ) J3w3 3 

\ MA4 4 — 5'(u;4 4 ) J4u?4 4 / 


Tfif,c — Fn,c 5(o,c c ) c — 1) 2, 3, 4, cl, c2, t 

o>c c = (<£ c , 4cos^ c , — 0 c sin<£ c ) T c = cl,c2 

wt = i3 -J- o,3 


( 1 —sin<f> 3 tan 0 3 —cos<j> 3 tan 6 3 \ / <}>t 

0 cosA<f> sinA(f> 1 j o>3 3 (2) 

0 —sinA 4 > coaA<t> J ya>3 3 (3) 

A 25 v\\ = -T n ,i 5'(u>l 1 )5'(i21*li) = T n , 1 5 2 (o,l!) JZl’li 

A 2 g o>2 2 = —Tnj S 2 (w22)R2*22 

A27 o>3 3 = — T^ )3 £ 2 (o>3 3 )(i23*3 3 — R3*4 3 ) 

A 3T o,3 3 = -Tiv.s S 2 (u>3 3 )R3*3 3 

A 47 w3 3 = -Tn .3 5 2 (u,3 3 ) R3*3 3 + T Nl t S(o>t t ) BT o>3 3 + r N)t 5T o,3 3 
A 4 $ o,4 4 = —Fn,4 *9 2 (o,4 4 ).R4*54 
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Table 2. (continued) 
(e) Interaction force. 


/-FCljvX 

/ Ajcljv 

0 

0 

0 ^ 


FC2 n 


0 

kcltf 

0 

0 


fcz n 


— A:cl tv 

—kc 2jv 

kc3w 

kcApf 

/TCI 

FC 4 N 

= H TC = 

0 

0 

—kcZiv 

—kc 4;v 

( TC2 

MCU 

^lli 

0 

0 

0 

TC3 

MC 2 2 


0 

£22 2 

0 

0 

Vrc4 

MC 3 3 


-*31, 

-£32 3 

£333 

£34 3 


M C4 4 / 


V 0 

0 

CO 

1 

-e44 4 v 



/ fccljv 

0 

0 

0 N 

0 

kc 

0 

0 

— &cl ^ 

— kc 2 pj 

*3jv 

ktpj 

0 

0 

-*3jv 

—ktpj 

aii 

0 

0 

0 

0 

£22 2 

0 

0 

-ai 3 

-a2 3 

a® 3 

(Zz 3 

0 

0 

-a®4 

-£4*4 / 


/TCI \ 

rc2 

FTx 
FTz / 


TC1,...!TC4 = cable tensions for Cl,. ..(74 


Ft = 

: kc3 TC3 

+ 

kc4 TC 4 

= i3 FTx + kt FTz 






ai 

= Rl*l 

X 

kcl 

II 

CO 

CO 

R3*3 

X 

kc3 

II 

X 

CO 

VJ> 

f3z jt 



= R3*3 

X 

kcl 

II 

CO 

T}1 

R3*5 

X 

kc3 

*Tirv 

X 

II 

R4*5 x 

i3 

^22 

= R2*2 

X 

kc2 

II 

CO 

w* 

R3*4 

X 

kc4 

£3z = 

-iZx jt 


M 

eo 

Vy» 

= R3*4 

X 

kc2 

II 

vv> 

R4*5 

X 

kc4 

II 

N 

R4*5 x 

kt 


i?13Ar = £1 A;cl at = / T'at.c! V^13 e i dt 

R24ff = £2 kclf] = fTff iC 2 y24 c2 dt 

R3*5* = / 7V >t ( VTt + A 7r f ) dt 

R35n — £3 kc3ff = R3*5 n + jT^p S(R3*3 3 ) u;3 3 dt 

R45fj = i 4 kc4 N = R3*5n + fTfj^ S(R3*4$) m3 3 dt 
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Table 2. (continued) 
Elastic cable tensions 


TCj = max{ 0, Kj(£j - £oj)} j = 1, ...,4 
Inelastic cable constraint force parameters 
s = -[A t D~ x A]" 1 A t D~ x fo 

Ai T D -, /o = kcl l (FOljv/mi - F03 N /m 3 ) + - &l 3 T J3~ x M03 3 

A 2 T D -1 /o = (F02 N /m 2 - F03 N /m 3 ) + ^22 2 T J2" 1 M02 2 - £32jj3 _1 M03 3 

A^D-'fo = t'3^ (F03*/m 3 - F04 N /m 4 ) + £3x 3 T J3" 1 M03 3 - £4x2J4 _1 M04 4 
A 4 t JD- , /o = (F03 N /m 3 - F04 N /m 4 ) + ^3z 3 T J3~ X M03 3 - £4zJj4 -1 M04 4 

Ai T 2?" 1 Ai =Mi3+^lli T Jl _1 eili +£3lfJ3- , £31 3 
A 2 T D~ X A 2 = /i 23 + £22 2 T J2 _1 £22 2 + ^32^J3-^32 3 
| A 3 T D~ X A 3 = /i 34 + £3x 3 T J3 -1 £3x 3 + £4xf J4 -1 £4x 4 

1 A 4 T D" 1 A 4 = /x 34 + £3z 3 T JS- i ^z 3 + {4*J’J4- 1 ^4z 4 

A i T D~'A 2 = Jfecl£ kc2 N /m 3 + ^31 3 T J3-^32 3 

I 

Ai T D~ x A 3 = —(kcljf iZ N /m 3 + £31 3 T J3~^3* 3 ) 

A 1 t D" 1 A 4 = — (fccljv ktN/m 3 + |31 3 t J3-^3z 3 ) 

A 2 t D~ x A 3 = —(kc2fj i3^/m 3 + £32 3 T J3~ *£3x 3 ) 

A 2 t D _1 A 4 = — (kc2pf ktpj/m 3 + £32 3 r J3 -1 £3z 3 ) 

A 3 t D- 1 A 4 = ^3x 3 t J3-^3z 3 + e4xfJ4-^4z 4 

f 

where 

/*13 = ( m i + m 3 )/m 1 m 3 
^23 = (”*2 + m 3 )/ m 2 m 3 
^34 = (ui 3 4- m 4 )/m 3 m 4 
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Table 2. (concluded) 

(f) Simulation equations. 


sf ^D-'ifo + fc) 


/sfl N \ 
s/2n 
sf 3w 
s/4n 
sm li 
am. 2 2 
am33 
\ am44 / 


/ (FOljy + jPCIjv)^ \ 

(F02jv + FC2N)/m2 
(F03 n + FCSn)/™* 
(F04jv + FC4 N )/m 4 
Jl-^MOh+MCh) 
J2" 1 (M02 2 + MC2 2 ) 
J3 -1 (M03 3 + MCZ%) 
\J4- 1 {M04 4 +MC4 4 )J 


ti = A 1 a/ 

Vl^ = a/liv 
u>li = amli 
u>2 2 = am2 2 
w3 3 = srrt&z 
u>4 4 = sm4 4 

F13 c i = T c i t N[sf3N ~ sflN ~ A 25 sm\\ - A 37 smSs] 
V"24 c2 = T c2) ]v[a/3jv ~ «/2j\r + A 2 g am2 2 — A 27 am3 3 ] 
VT t — T t) N[sf4N — s/Zn ~ A 4 8 am4 4 ] — BT am3 3 
Note: 

/ .<* T \ 

T C ,N = JCN T ) c — cl,c2,< 

\kc N T ) 
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Secondly, the constraint force on the configuration, /c, is given in table 2e. In this example, equa- 
tion (11) or equations (18, 19) can be used for both elastic and inelastic cables, because c = m. Equa- 
tion (11) gives fc in terms of the four cable tensions and the coefficient matrix, H, given in table 3e. 
Here, H is also a basis of the constraint force space. A second basis is given formally by the rows of A~ x 
corresponding to the cable stretch rates (rows 6, 9, 10, 12) and this is noted as A in the table 2e. The 
corresponding constraint force parameters can be identified (by equating the two expressions for fc) as 
the tether cable tensions, TC\,TC2, and the ? t components of the constraint force, FT acting on the 
triangle at R5 ; that is, the F t components of the sum of bridle cable forces. Formulas for calculating cable 
line segments from u are included in the appendix, and these suffice to determine H and A , as well as 
{TCj, FTt }, for elastic cables. For inelastic cables, the constraint force parameters corresponding to A 
or H are given by equation (19). Expressions for the scalar elements of equation (19) corresponding to A 
are listed in the table. Similar expressions corresponding to H can be routinely generated, and the required 
inverse of the 4 x 4 matrix, [ A r D~ l A], can be left to ordinary numerical methods without difficulty. 

Finally, the simulation equations are assembled in table 2f. First the total specific force on the config- 
uration, sf, due to external forces, fo (which define the free body motion), and the interaction due to the 
suspensions, fc , are given. The vector elements of sf are denoted sfl , sm4 . The simulation equations 
are given by expanding the relation u = A~ x sf, representing equations (16) or (20). For inelastic cables, 
the cable stretching DOFs, \ , (the components of V 13 along kcl , V 24 along kc2 , and V T along i3 , kt ) 
are all theoretically zero. 

Table 2 provides programmable simulation equations for a rigid body model of the dual lift system 
with spreader bar for elastic or inelastic cables. The expressions for fo, fc, u can be rearranged further 
to reduce computational requirements. Repeated matrix-vector combinations are apparent in these expres- 
sions and are readily eliminated. Cross products occur frequently in the vector elements of Au, A , and an 
efficient generic routine for this is useful. Programming in terms of the matrix and vector objects of rigid 
body mechanics used in table 2 is possible with the language ADA ([35]), together with the matrix algebra 
package of [36]. 

The other dual lift systems, 5 and 7, of figure 1 are both three-body subsystems of system 6 obtained 
by deleting the load and bridle cables. System 5 can be represented by the same generalized velocity 
coordinates and equations as system 6 after deleting the six load-triangle coordinates and corresponding 
load forces and moments. System 7 is a specialization of system 5 with both cables attached at the same 
point on the load. Thus, all three dual lift systems are readily integrated in a single simulation. 

5. CONCLUSIONS 

General simulation equations for slung load systems, derived from D’Alembert’s principle, account for 
any suspension geometry including controllable geometry, cable elasticity, and choice of generalized co- 
ordinates. These equations generalize previous results for single helicopter systems and provide a new 
formulation for the case of inelastic suspensions, which is computationally more efficient than the conven- 
tional formulation, unifies the elastic and inelastic cable formulations, and is readily applied to the more 
complex dual lift and multilift systems. 

Application of the general equations to the derivation of simulation equations for specific systems is 
demonstrated in the case of dual lift systems, and this previously difficult problem is seen to be tractable 
with these methods. 
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APPENDIX: SKEW-SYMMETRIC MATRICES AND 
COORDINATE TRANSFORMATIONS 

The working equations of this report are given entirely in terms of vector functions from three-dimensional 
rigid body dynamics. These are expressed in a programmable form with the aid of skew-symmetric ma- 
trices and coordinate transformations. The general skew- symmetric matrix S( x, y, z) is defined from the 
scalar triplet (x, y, z) as given in table Al. This allows scalar representation of the vector cross products 
(VI x V2 ) as shown in the table. The algebra of skew-symmetric matrices is consistent with correspond- 
ing relations from vector algebra, such as sums and reversals of cross products and triple product relations. 
Cross products pervade the dynamics in coriolis velocities and accelerations and centrifugal accelerations. 
Corresponding scalar expressions are noted in the table. Cross products representing moments due to ca- 
bles occur in the equations of motion. A standardized formulation of this moment used in the present paper 
is given in the table. If cable Cj applies tension TCj in the direction kcj at point Rj on body £,, then its 
moment about the c.g. of B,- is given by 

Mij = Ri*j x kcj TCj = £ij TCj 

where the symbol £ij is reserved for the moment action of cable Cj on body £, per unit tension. 

Euler angle transformations appear in all dynamic terms, and define the transformation of vectors given 
by their coordinates in inertial space, to their coordinates in body or cable axes in terms of a sequence of 
single axis rotations. The primitive transformations for rotations about a single axis i ,j , k of a right-handed 
orthogonal reference frame are denoted E\ ( ct) , Ez ( a) , Ej ( cr) , and are given in table A2. Transformations 
to body axes T^n are given by the usual yaw, pitch, and roll rotation sequence of aeronautical usage, and 
are listed in table A2. For this report, cable axes are defined from inertial pitch and roll angles which locate 
the cable direction. This is a truncated Euler sequence and the corresponding transformation is given by 
nulling the yaw angle in Tj,# • In a simulation, Euler angles can be obtained from the kinematic angular 
rate relation ctb(ubb) given in the table. Transformation rates occur throughout the inertia reaction term, 
Au, which represents coriolis and centrifugal accelerations. The coordinates of an arbitrary vector, V , in 
frames jFb, ?n are related by 

V N = T NJb V b 

The derivative 

Vn - Tff,b Vb + TV , b Vb 

is the scalar counterpart of the usual vector coriolis equation relating time derivatives relative to two dif- 
ferent coordinate frames. The second term can be identified as the effect of coordinate frame rotation, so 
that 

T a , b Vb = (wbxV) 4 = T a>b S(u)b b )Vb 

and this is valid for arbitrary Vj. 
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Table Al. Skew-symmetric matrices and cross products 


( 0 -z y\ 

S(x,y,z) = 2 0 -x 

\ -y x 0 / 

(VI x V2 ) 0 = 5(Vl 0 )V2 a = — (V2 x VI ) 0 = -S(V2 0 )Vl a 


Coriolis and centrifugal terms 


(to x R) 0 = S(to a )R a = -S(R a )to a 
(w x V) 0 = S(u a )V a = —S(V a )to a 

(to xu/ xR) a = S 2 (to a )R a = -S(to a )S(R a )to a = w*Rw a - to 2 R a 


Moment of cable Cj on body Bi about Ri * 


Ti = body axes fcsrBi 

Miji = (Ri’j x rC;kcj), = S(Ri*j i )kcj i TCj = iijCTCj 
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Table A2. Coordinate transformations 


Single axis transformations 


/I 0 0 \ / cos ct 1 

E\ (ct) = I 0 cos a sin a J £^(ct) = I 0 1 

\ 0 — sin a cos a / V sin a 0 


- sin a \ / cos a sin a 0 \ 

0 ] Ei(cr) = I — sin ct cos cr 0 1 

cos a / \0 0 1 / 


Euler angle transformation. Th v 

J? N = {iN , jN , kN } = inertial axes 
JF, = {ib,jb,kb} = body axes 
a 6 = (0s, 6b, \pb) T = Euler angle triplet 





Tb'Ni&b) = Ei((f>b) E2(6b) Ez(ij)b) - 

cos 0s cos 6b sin ip b cos 6 b - sin 6 b 

sin <pb cos 0s sin 6b — cos 0s sin 0 & sin 0s sin 0 6 sin 0s + cos 0s cos 0s sin 0s cos 0s 
cos <pb cos 0s sin 0s + sin 0s sin 0 b cos 0s sin 0s sin 0s — sin 0s cos 0s cos 0s cos 0s . 


Angular velocity of Ek relative to jEu 

/I 0 - sin 6 b \ / 0s \ / 06 \ 

u)bb = PF(a£>)dfc = I 0 cos 0 4 sin 0 4 cos 6 b ) ( 06 ] = ( jb' h , kN b ) I 06 1 

Vo - sin 0s cos (f>b cos 0 6 / V 0s / V 0s / 


/I 

sin (f>b tan 0s 

cos 06 tan 0s \ 

/ ib ' • wb / cos \ 

&& = VF l (ab)u)bb = ( 0 

COS <f>b 

- sin (f> b wis = j 

jb' • wb 

Vo 

sin cos 0 6 

cos 0(,/ cos 6 b J 

V kb ' o wb / cos 6 b J 


Transformation rates 


TV, 6 = TV,6 S(wbb) 
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